/***************************************************************************

 Ulam Spiral Generator
 Copyright (C) 2012 Nick Lowery (nick.a.lowery@gmail.com)

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 www.gnu.org/copyleft/gpl.html

 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define ULAMSPIRAL_TEXT_BYTEWIDTH 8
#define ULAMSPIRAL_TEXT_BYTEWIDTH_STRING "8"

int numDivisors(unsigned int n);
int isPrime(unsigned int n);
int isPerfectNExp(unsigned int n, unsigned int exp);
void createGrid(void);
void writeGridBMP(void);
void cleanUp(void);

unsigned int ** g_numberGrid;
unsigned int g_gridSize;
unsigned int g_gridSizeSq;
enum {
	RENDERINGSTYLE_PRIMES, RENDERINGSTYLE_COMPOSITENESS, RENDERINGSTYLE_PERFECT,
} g_renderingStyle;
unsigned int g_exponent;

int main(void) {
	int okay;

	// get spiral size
	okay = 0;
	do {
		printf("\nSize of Ulam Spiral to create:\n");
		scanf("%u", &g_gridSize);

		// make sure the spiral size is valid
		if (g_gridSize == 0 || g_gridSize > 65535 || g_gridSize % 4 != 0) {
			printf(
					"%u is an invalid size. Must be a positive integer, less than 65,536, and a multiple of four.\n",
					g_gridSize);
		} else {
			g_gridSizeSq = g_gridSize * g_gridSize;
			okay = 1;
		}
	} while (!okay);

	// get spiral style
	printf("\nSelect a rendering style:\n"
			"1. Highlight prime numbers only\n"
			"2. Highlight every number based on compositeness\n"
			"3. Highlight perfect N-exponents only\n");
	scanf("%u", &g_renderingStyle);
	g_renderingStyle--;

	// get exponent if user wants perfect n-exponents
	if (g_renderingStyle == RENDERINGSTYLE_PERFECT) {
		okay = 0;
		do {
			printf(
					"\nEnter exponent (perfect squares=2, perfect cubes=3, etc):\n");
			scanf("%u", &g_exponent);

			if (g_exponent <= 1)
				printf(
						"%u is an invalid exponent. Must be greater than one.\n", g_exponent);
			else
				okay = 1;
		} while (!okay);
	}

	printf("\nWorking...\n");

	// create number grid and ulam spiral
	createGrid();

	// output spiral
	writeGridBMP();

	printf("Done!\n");

	// clean up memory and exit
	cleanUp();
	return 0;
}

void createGrid(void) {
	unsigned int val;
	unsigned int x, y;
	int dx, dy;
	unsigned int amt; // length of spiral's current side
	unsigned int mvd; // amount move along spiral's current side
	unsigned int c; // current number of direction switches

	// create grid in memory

	g_numberGrid = (unsigned int **) malloc(
			g_gridSize * sizeof(unsigned int *)); // 4 bytes per pointer
	for (x = 0; x < g_gridSize; x++) {
		g_numberGrid[x] = (unsigned int *) malloc(g_gridSize * sizeof(unsigned int)); // 4 bytes per number value
	}

	// fill in grid numbers

	// find starting point in grid
	if (g_gridSize % 2 == 0)
		x = y = (unsigned int) ((g_gridSize - 1) / 2.0f);
	else
		x = y = (unsigned int) (g_gridSize / 2.0f);

	dx = 1;
	dy = 0;
	val = 1;
	amt = 1;
	c = 0;
	while (val <= g_gridSizeSq) {
		mvd = 0;

		// move in the current direction, filling in values along the way
		while (mvd < amt && val <= g_gridSizeSq) {
			g_numberGrid[x][y] = val;
			x = (unsigned int) ((int) x + dx);
			y = (unsigned int) ((int) y + dy);
			mvd++;
			val++;
		}

		// every two switches in direction, the sprial's size increases by 1
		c++;
		if (c == 2) {
			c = 0;
			amt++;
		}

		// cycle directions to move around the spiral in
		if (dx == 1) {
			dx = 0;
			dy = 1;
		} else if (dy == 1) {
			dy = 0;
			dx = -1;
		} else if (dx == -1) {
			dx = 0;
			dy = -1;
		} else if (dy == -1) {
			dy = 0;
			dx = 1;
		}
	}
}

void writeGridBMP(void) {
	// base code from http://stackoverflow.com/questions/2654480/writing-bmp-image-in-pure-c-c-without-other-libraries

	int x, y;
	FILE * out;
	int w = (int) g_gridSize;
	int h = w;
	unsigned char * img = NULL;
	const int filesize = 54 + 3 * w * h;
	unsigned char r, g, b;
	float intensity;
	unsigned char bmpfileheader[14] = { 'B', 'M', 0, 0, 0, 0, 0, 0, 0, 0, 54, 0,
			0, 0 };
	unsigned char bmpinfoheader[40] = { 40, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
			0, 24, 0 };
	unsigned char bmppad[3] = { 0, 0, 0 };
	char filename[64] = { 0 };
	char gridsize[8] = { 0 };

	img = (unsigned char *) malloc(filesize - 54); // BMP headers are 54 bytes

	// only show primes
	if (g_renderingStyle == RENDERINGSTYLE_PRIMES) {
		for (x = 0; x < w; x++) {
			for (y = 0; y < h; y++) {
				if (isPrime(g_numberGrid[x][y])) {
					r = g = b = 255;
				} else {
					r = g = b = 0;
				}

				img[(x + y * w) * 3 + 2] = r;
				img[(x + y * w) * 3 + 1] = g;
				img[(x + y * w) * 3 + 0] = b;
			}
		}
		// show all integers where color brightness indicates number compositeness
	} else if (g_renderingStyle == RENDERINGSTYLE_COMPOSITENESS) {
		for (x = 0; x < w; x++) {
			for (y = 0; y < h; y++) {
				intensity = numDivisors(g_numberGrid[x][y]) * 4;

				img[(x + y * w) * 3 + 2] = (unsigned char) (intensity);
				img[(x + y * w) * 3 + 1] = (unsigned char) (intensity);
				img[(x + y * w) * 3 + 0] = (unsigned char) (intensity);
			}
		}
		// highlight perfect n-exponents
	} else if (g_renderingStyle == RENDERINGSTYLE_PERFECT) {
		for (x = 0; x < w; x++) {
			for (y = 0; y < h; y++) {
				if (isPerfectNExp(g_numberGrid[x][y], g_exponent)) {
					r = g = b = 255;
				} else {
					r = g = b = 0;
				}

				img[(x + y * w) * 3 + 2] = r;
				img[(x + y * w) * 3 + 1] = g;
				img[(x + y * w) * 3 + 0] = b;
			}
		}
	}

	// BMP is little-endian, so reverse byte order
	bmpfileheader[2] = (unsigned char) (filesize);
	bmpfileheader[3] = (unsigned char) (filesize >> 8);
	bmpfileheader[4] = (unsigned char) (filesize >> 16);
	bmpfileheader[5] = (unsigned char) (filesize >> 24);
	bmpinfoheader[4] = (unsigned char) (w);
	bmpinfoheader[5] = (unsigned char) (w >> 8);
	bmpinfoheader[6] = (unsigned char) (w >> 16);
	bmpinfoheader[7] = (unsigned char) (w >> 24);
	bmpinfoheader[8] = (unsigned char) (h);
	bmpinfoheader[9] = (unsigned char) (h >> 8);
	bmpinfoheader[10] = (unsigned char) (h >> 16);
	bmpinfoheader[11] = (unsigned char) (h >> 24);

	// generate filename

	sprintf(gridsize, "%u", g_gridSize);
	strcat(filename, "ulamspiral_");
	strcat(filename, gridsize);
	if (g_renderingStyle == RENDERINGSTYLE_PRIMES)
		strcat(filename, "_p");
	else if (g_renderingStyle == RENDERINGSTYLE_COMPOSITENESS)
		strcat(filename, "_c");
	else if (g_renderingStyle == RENDERINGSTYLE_PERFECT)
		strcat(filename, "_pe");
	strcat(filename, ".bmp");

	// create file

	out = fopen(filename, "w");
	if (out == NULL) {
		printf("Failed to save spiral to BMP.\n");
		free(img);
		return;
	}

	// write BMP to file

	fwrite(bmpfileheader, 1, 14, out);
	fwrite(bmpinfoheader, 1, 40, out);
	for (y = 0; y < h; y++) {
		fwrite(img + (w * (h - y - 1) * 3), 3, w, out);
		fwrite(bmppad, 1, (4 - (w * 3) % 4) % 4, out);
	}
	fclose(out);
	free(img);
}

void cleanUp(void) {
	unsigned int x;

	if (g_numberGrid != NULL) {
		for (x = 0; x < g_gridSize; x++) {
			if (g_numberGrid[x] != NULL)
				free(g_numberGrid[x]);
		}

		free(g_numberGrid);
	}
}

int isPrime(unsigned int n) {
	unsigned int divisor;
	const unsigned int bound = (unsigned int) (sqrt((float) n) - 1);

	if (n <= 1 || n % 2 == 0)
		return 0;

	for (divisor = 3; divisor <= bound; divisor += 2) {
		if (n % divisor == 0)
			return 0;
	}

	return 1;
}

int numDivisors(unsigned int n) {
	unsigned int divisors = 0;
	unsigned int divisor;;
	const unsigned int bound = (unsigned int) (sqrt((float) n));

	for (divisor = 3; divisor <= bound; divisor++) {
		if (n % divisor == 0)
			divisors++;
	}

	return divisors;
}

int isPerfectNExp(unsigned int n, unsigned int exp) {
	const float root = (float) pow((double) n, 1.0 / (double) exp);
	const float rootfloor = floor((double) root);

	return ((root - rootfloor) == 0.0f);
}
